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Abstract 

This paper presents a brief overview of some recent ad- 
vances in numerical radiative transfer, which may help the 
molecular astrophysics community to achieve new break- 
throughs in the interpretation of spectro-(polarimetric) 
observations. 
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1. Introduction 

The development of novel Radiative Transfer (RT) meth- 
ods often leads to important breakthroughs in astrophysi- 
cal plasma spectroscopy because they allow the investiga- 
tion of problems that could not be properly tackled using 
the methods previously available. This RT topic is also of 
great interest for the "Promise of FIRST" , since a rigor- 
ous interpretation of the observations will require to carry 
out detailed confrontations with the results from NLTE 
RT simulations in one-, two-, and three-dimensional ge- 
ometries. 

The efficient solution of NLTE multilevel RT problems 
requires the combination of a highly convergent iterative 
scheme with a very fast formal solver of the RT equation. 
This applies to the case of unpolarized radiation in atomic 
lines, to the promising topic of the generation and transfer 
of polarized radiation in magnetized plasmas and to RT 
in molecular lines. 

The "dream" of numerical RT is to develop iterative 
methods where everything goes as simply as with the well- 
known A-iteration scheme, but for which the convergence 
rate is extremely high. In this contribution we present 
an overview of some iterative methods and formal solvers 
we have developed for RT applications. Our RT methods 
are based on Gauss-Seidel iteration and on the non-linear 
multigrid method. These new RT developments are of in- 
terest because they allow the solution of a given RT prob- 
lem with an order-of-magnitude of improvement in the to- 
tal computational work with respect to the popular ALI 
method (on which most present NLTE codes are based 
on) . Our RT methods have been succesfuUy applied to as- 
trophysical problems of unpolarized radiation in atomic 
lines (in ID, 2D and 3D with multilevel atoms) and also 



to the transfer of polarized radiation in magnetized plas- 
mas including anisotropic pumping (Trujillo Bueno, 1999; 
2001), which may be of interest for modelling polarization 
phenomena in MASERS. The case of RT in molecular lines 
is presented in an extra contribution at this conference by 
Asensio Ramos, Trujillo Bueno and Cernicharo (2001). 

2. RT METHODS BASED ON GAUSS-SeIDEL ITERATION 

The essential ideas behind the iterative schemes on which 
our NLTE multilevel transfer codes are based on can be 
easily understood by considering the "simplest" NLTE 
problem: the coherent scattering case with a source func- 
tion given by 



(l-e)J + eB, 



(1) 



with e the NLTE parameter, J the mean intensity and 
B the Planck function. The mean intensity at point "z" 
is the angular average of incoming ( "in" ) and outgoing 
("owt") contributions. For instance, for a one-point angu- 
lar quadrature 



Ji = Ji 
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(c+r*)- 



(2) 



The well-known A— iteration scheme is to do the fol- 
lowing in order to obtain the "new" estimate of the source 
function at each spatial grid-point "i" : 



.(l-e)jr + eBi, (3) 

where J°^'^ is the mean intensity at the grid-point "i" cal- 
culated using the previous known values of the source 
function (i.e. using Sf'^). 

For a given spatial grid of NP points the formal so- 
lution of the transfer equation can be symbolically repre- 
sented as 



Ir 



An [S] 



(4) 



where To gives the transmitted specific intensity due to 
the incident radiation at the boundary and A^ is a NP x NP 
operator whose elements depend on the optical distances 
between the grid-points. Thus, the mean intensity at the 
grid-point "j" would be: 



Ai,iS? + ... 
^Ai^i+iS^_^j 



Ai,i-iSf_]^ + AijSj + 
... + Aj^NpS^p + Ti. 



(5) 
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In this last expression A^- = E(A*"+A™*)/2 (with the 
sum appUed to all the directions of the numerical angular 
quadrature) and a, b and c are simply symbols that we use 
as a notational trick to indicate below whether we choose 
the "old" or the "new" values of the source function. Thus, 
for instance, the A-iteration method consists in calculating 



Told 



as 



Ji choosing a = b = c = old, which gives Ji 
indicated in Eq. (3). 

The Jacobi method, known in the RT literature as the 
OAB method (from Olson, Auer and Buchler, 1986), and 
on which most NLTE codes are based on, is found by 
choosing a = c = old, but b = new, which yields 

Ji = J?'^ + Aii (Sr^ - S?'<i) = J?''^ + Aii 6Si (6) 

In fact, using this expression instead of Jf'^ in Eq. (3) 
one finds that the resulting Jacobi iterative scheme is 



with the correction 



<5Si 



[(1 - e) J°id + eBj - SP"] 
[l-(l-^)Aii] ■ 



(7) 



(8) 



where An is the diagonal element of the A-operator as- 
sociated to the spatial grid-point "i" . Note that the cor- 
rection corresponding to the slowly convergent A-iteration 
method is given by Eq. (8), but with An = 0. 

A superior type of iterative schemes are the Gauss- 
Seidel (GS) based methods of Trujillo Bueno and Fabi- 
ani Bendicho (1995), which can also be suitably general- 
ized to the polarization transfer case (cf. Trujillo Bueno 
& Landi Degl'Innocenti, 1996; Trujillo Bueno and Manso 
Sainz, 1999). This type of iterative schemes are obtained 
by choosing c = old and a = b = new. This yields 



Ji = j: 



old and new 



+ An 6Si 



(9) 



where jpi^andnew jg ^jjg mgau intensity calculated using 

the "new" values of the source function at grid-points 
l,2,...,z— 1 and the "old" values at points z, z-t- 1, i -|- 2, 
NP. The iterative correction is given by 



5S?^ = 



[(1-c ).];'•' 



laudi 



[l-(l-e)Aii] 



(10) 



It is important to clarify the meaning of this last equa- 
tion: 

1) First, at point i = 1 (which we can freely be as- 
signed to any of the two boundaries of the medium under 
consideration) use "old" source function values to calcu- 
late Ji via a formal solution. Apply Eqs. (10) and (7) to 
calculate Sf ^. 

2) Go to the next point i = 2 and use S"*'^ and the 
"old" source-function values S"'*^ at points j = 2, 3, ...,NP 
to get J2 via a formal solution. Apply Eqs. (10) and (7) 
to calculate 

3) Go to the next spatial point k and use the previously 

obtained "new" source fimction values at j = 1, 2, fc — 1, 
but the still "old" ones at j = k,k + l, ...,NP to get J^ via 



a formal solution and Sg**"" as dictated by Eqs. (10) and 
(7). 

4) Go to the next point k + 1 and repeat what has 

just been indicated in the previous point until arriving 
to the other "boundary point" . Having reached this stage 
iniciate again the same process, but choosing now as "first 
point i = 1" the above-mentioned "boundary point" . 

The result of what we have just indicated is a pure GS 
iteration. Actually, after an incoming and outgoing pass 
we get two GS iterations! A SOR iteration is obtained by 
doing the corrections as follows: 

5SfO^ = a;(5Sps, (11) 

where w is a parameter with an optimal value between 1 
and 2 that can be found easily (see Trujillo Bueno and 
Fabiani Bendicho, 1995). 

Figure 1 shows an example of the convergence rate of 
all these iterative methods applied to a NLTE polarization 
transfer problem in a stellar model atmosphere permeated 
by a constant magnetic field that produces a Zeeman split- 
ting of 3 Doppler widths. We point out that the computing 
time per iteration is similar for all these methods and that 
matrix inversions are not performed. Thus, the important 
point to keep in mind is that our implementation of the 
GS method is 4 times faster than Jacobi (i.e. than the ALI 
method on which most NLTE codes are based on), while 
our SOR method for radiative transfer applications is 10 
times faster. 
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Figure 1. The variation of the maximum relative change 
versus the iteration number for several types of iterative 
methods applied to the NLTE Zeeman line transfer prob- 
lem discussed by Trujillo Bueno and Landi Degl'Innocenti 
(1996). The NLTE parameter e = 10"''. Dotted line: 
the A-iteration method. Dashed line: the Jacohi-based ALI 
method. Solid line: the GS-based method. Dashed-dotted 
line: the SOR method. 



For pedagogical reasons we have chosen here a NLTE 

linear problem in order to explain in simple terms our 
GS-based iterative schemes. The generalization to the full 
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non-linear multilevel problem can be carried out as in- 
dicated in the Appendix of the paper by Trujillo Bueno 
(1999). The critical point is always to remember that the 
approximations one introduces for achieving the required 
linearity of the statistical equilibrium equations at each it- 
erative step should treat adequately the coupling between 
transitions and the non-locality of the problem (see Socas- 
Navarro & Trujillo Bueno, 1997). 

3. The non-linear multigrid RT method 

Our GS and SOR radiative transfer methods are based, 

like the Jacobi-likc ALT method, on the idea of operator 
splitting. Therefore, all of them are characterized by a 
convergent rate which decreases as the resolution of the 
spatial grid is increased. As a result, if NP is the number 
of grid-points in a computational box of fixed dimensions, 
the computing time or computational work (W) required 
by the three previous iterative methods to yield the self- 
consistent atomic (or molecular) level populations scales 
with NP as follows (cf. Trujillo Bueno & Fabiani Bendicho, 
1995): 

- Jacobi-based ALT method yV^-NP^ 

- Our GS-based RT method ^ >V~NPV4 

- Our SOR RT method ^ W~NP\/NP 

Is there any suitable RT multilevel method character- 
ized by yy--NP? This would be of great interest for 3D RT 
with multilevel atoms where NP~ 10^. The answer is af- 
firmative. This has been worked out by Fabiani Bendicho, 
Trujillo Bueno and Auer (1997) who considered the appli- 
cation of the non-linear multigrid method (see Hackbush, 
1985) to multilevel RT. 

The iterative scheme of the non-linear multigrid method 
is composed of two parts: a smoothing one where a small 
number of GS-based iterations on the desired finest grid 
are used to get rid of the high-frequency spatial compo- 
nents of the error in the current estimate, and a correction 
obtained from the solution of an error equation in a coarser 
grid. With our non-linear multigrid code the total com- 
putational work scales simply as NP, although it must be 
said that the computing time per iteration is about 4 times 
larger than that required by the A— iteration method. 

In order to compare the convergence rate of all these 
iterative methods we present in Fig. 2 an estimate of 
the maximum eigenvalues (p) of the corresponding iter- 
ation operator, which controls the convergence properties 
of such iterative schemes. The knowledge of this maxi- 
mum eigenvalue (p) is usehil because errors decrease as 
p***", where "ztr" is the iterative step. As it can be noted in 
Fig. 2 the convergence rate of both, the MALI and MUGA 
schemes decreases when the spatial resolution of the grid is 
improved, while the maximum eigenvalue of our non-linear 
multigrid method is always very small (p 0.1) and in- 
sensitive to the grid-size. A maximum eigenvalue p = 0.1 
means that the error decreases by one order of magnitude 
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Figure 2. Variation with the grid-spacing Az of the maxi- 
mum eigenvalue of the iteration operator corresponding to 
several multilevel iterative schemes. The MG symbol refers 
to our non-linear multigrid code, while MU GA to our mul- 
tilevel GS-based code. MALI refers to the Jacobi-based 
multilevel ALI method of Rybicki & Hummer, (1991). 

each time we perform an iteration! This explains that, typ- 
ically, two multigrid iterations are sufficient to reach the 
self-consistent solution for the atomic level populations. 

4. Formal solvers for RT applications 

The formal solution routines of our NLTE codes (for unpo- 
larized or polarized radiation and for atomic or molecular 
species) are based on improvements and generalizations 
of the short- characteristics (SC) technique introduced by 
Kunasz & Auer (1988). Let us recall it briefly indicating 
also our generalization to 3D radiative transfer and to the 
case of polarized radiation. 

The scalar RT equation for the specific intensity is 

^=X.(S. -I.), (12) 

where s is the geometric distance along the ray propagat- 
ing in a certain direction in a 3D medium, Xv is the total 
opacity and Si^ the source function. 

Point O is the grid-point of interest at which one wishes 
to calculate the specific intensity Iq, for a given frequency 
[v) and a direction (fJ). Point M is the the intersection 
point with the grid-plane that one finds when moving 
along — n. At this upwind point the specific intensity Im 
(for the same frequency and angle) is known from pre- 
vious steps. In a similar way, point P is the intersection 
point with the grid-plane that one encounters when mov- 
ing along Jl. We also introduce the optical depths along 
the ray between points M and O (Atm) and between 
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points O and P (Arp). From the formal solution of the 
previous transfer equation one finds that 

/■Atm 

Io=lMe-^^+/ S(i)e-(^™-*)rfi, (13) 
Jo 

with the optical depth variable measured from M to O. 

The integral of this equation can be solved analytically 
by integrating along the short- characteristics MO assum- 
ing that the source function S{t) varies parabolically along 
M,0 and P. The result reads: 

lo = iMe-^^ + *mSm + *oSo + *pSp, (14) 

where (with X either M, O or P) are given in terms of 
the quantities Atm and Arp that we evaluate numerically 
by assuming that ln(x) varies linearly with the geometrical 
depth, X being the opacity. 

If the interest lies in the generation and transfer of po- 
larized radiation in magnetized astrophysical plasmas (cf. 
Trujillo Bueno and Landi Degl'Innocenti, 1996; Trujillo 
Bueno, 1999; 2001) the situation is a bit more complicated 
because, instead of having to solve the previous RT equa- 
tion for the specific intensity, one has to solve, in general, a 
vectorial transfer equation for the four Stokes parameters. 
For example, for the standard case of polarization induced 
by the Zeeman effect, the Stokes- vector at the grid-point 
O is 

/■Atm 

lo = 0(0, Atm)Im + / 0{t, Atm) S(i) dt, (15) 
Jo 

where 0{t,ATM) is the evolution operator (i.e. the 4x4 

Mueller matrix of the atmospheric slab between t and 
Atm)- In general, this evolution operator does not have 
an easy analytical expression and the integral of the previ- 
ous equation cannot be solved analytically. However, if the 
4x4 absorption matrix conmutes between depth points M 
and O (e.g. because one assumes the absorption matrix to 
be constant between M and O and equal to its true value 
at the middle point) the evolution operator reduces then 
to an expression given by the exponential of the absorp- 
tion matrix. The integral of Eq. (15) can then be solved 
analytically provided that the source function vector S is 
assumed to vary parabolically along M, O and P. Our for- 
mal solution method for Zeeman line transfer is precisely 
based on this idea and the results reads: 

lo = 0(0, ArM)lM + *mSm + *oSo + *pSp, (16) 

where (with X either M, O or P) are 4x4 matrices 
which are given in terms of Atm and Arp, in terms of 
the inverse of the absorption matrix and in terms of the 
analytical expression given by Landi Degl'Innocenti and 
Landi Degl'Innocenti (1985) for the evolution operator for 
the case in which the absorption matrix is assumed to be 
constant between the spatial grid points. An alternative 
formal solver of the Stokes- vector transfer equation suit- 
able for NLTE applications is the one given by Eqs. (1-4) 
of Socas-Navarro, Trujillo Bueno & Ruiz Cobo (2000). 



The application of these formal solution methods in 
ID is straightforward. The generalization to 2D geome- 
tries with horizontal periodic boundary conditions is sub- 
stantially more complicated. A suitable strategy has been 
described by Auer, Fabiani Bendicho and Trujillo Bueno 
(1994). 

The main changes when going to 3D imposing hori- 
zontal periodic boundary conditions lie in the interpola- 
tion. We have assumed that Im is known but, in most 
cases, the M-point (like the point P) will not be a grid- 
point of the chosen 3D spatial grid. The intensity at this 
M-point has to be calculated by interpolating from the 
available information at the nine surrounding grid-points, 
as we must also do for obtaining the opacities and source 
functions at M and P. Parabolic interpolation can how- 
ever generate spurious negative intensities if the spatial 
variation of the physical quantities is not well resolved 
by the spatial grid. This happens, for instance, if one 
tries to simulate the propagation of a beam in vacuum 
using a three dimensional grid. To avoid these problems 
we have improved the ID monotonic interpolation strat- 
egy of Auer and Paletou (1994), and generalized it to the 
two-dimensional parabolic interpolation that is required 
for 3D RT calculations with multilevel atomic models (see 
Fabiani Bendicho & TrujiUo Bueno, 1999). 

5. Conclusions 

The RT methods presented here are especially attractive 
because of their direct applicability to a variety of compli- 
cated RT problems of astrophysical interest. We empha- 
size that their convergence rate are extremely high, that 
they do not require the construction and the inversion of 
any large matrix and that the computing time per itera- 
tion is very small. 
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